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Abstract. In this article, we review the use of numerical techniques to obtain solutions for 
the quantum Hamiltonian constraint in loop quantum cosmology (LQC). First, we summa- 
rize the basic features of LQC, and describe features of the constraint equations to solve - 
generically, these are difference (rather than differential) equations. Important issues such 
as differing quantization methods, stability of the solutions, the semi-classical limit, and the 
relevance of lattice refinement in the difference equations are discussed. Finally, the cosmo- 
logical models already considered in the literature are listed, along with typical features in 
these models and open issues. 
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1 Introduction 

One of the main open questions for any feasible theory of quantum gravity is that of resolution 
of classical singularities. Due to the lack of experimental verification in this context currently, 
the resolution of classical singularities as well as the recovery of the appropriate semi-classical 
behavior are the key check that any acceptable theory of quantum gravity must overcome. 
Therefore, testing these two issues has been the main goal of all the numerical work performed 
in the symmetry reduced context of loop quantum cosmology (LQC). 

The geometrodynamical quantization of general relativity [48] is assumed to be the correct 
semi-classical limit of any quantum theory of gravity. In this context, the metric (or in the 
symmetry reduced cosmological context, the scale factor) is understood as the fundamental 
variable. Nevertheless, this quantization can not be fundamental because the quantum states 
completely follow the classical trajectories of the universe in the phase space. More precisely, 
if a state representing the current classical universe is evolved backwards making use of the 
Wheeler-DeWitt (WdW) equation, it ends up in the big bang singularity. As we will explain 
below, the situation in LQC is completely different and, in all analyzed models, the resolution 
of the corresponding singularity has been reported. 

This article is organized as follows: In Section 2, we provide a brief introduction to LQC with 
emphasis on the different quantization schemes used. Section 3 discusses (numerical analysis 

*This paper is a contribution to the Special Issue "Loop Quantum Gravity and Cosmology" . The full collection 
is available at http://www.emis.de/journals/SIGMA/LQGC.html 
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inspired) stability properties of the equations that arise in LQC models. In Section 4, we survey 
a number of isotropic and non-isotropic models. We end this article with a brief summary in 
Section 5. 

2 Essentials of loop quantum cosmology 

2.1 Phase space variables and quantum operators 

In this section, we will develop the general setting of LQC models - in particular, the origin 
of the Hamiltonian constraint difference equation. This will allow us to make some generic 
comments about solving these constraints using numerical techniques. In LQG, an oriented 
spatial manifold M is chosen, with a fiducial set of spatial triad vectors mapping from the 
internal vector space (with a positive definite metric qij) to the tangent space (with indices 
a,b, . . .) of the manifold M^. The dual of the triads are the 1-forms cj^, which give a positive 
definite metric Qab = Qij'^a'^l on M. For the Bianchi models, these 1-forms satisfy the Maurer- 
Cartan relation duj"^ = {l/2)Cjf^uj^ Auj^, where Cj^ are the structure constants of the Lie algebra 
for the particular Bianchi system. 

For the following, we consider specifically the case of a type I Bianchi model [9], but similar 
results hold for other systems as well. In particular, we use a diagonal form of Bianchi I, so that 
the physical metric is given in terms of three scale factors aj(T) and the lapse M by 

ds'^ = -M'^dr'^ + al{T)dxl + al{T)dxl + al{T)dxl. 

As with all Bianchi cosmologies, the spatial slices are non-compact and all physical quantities 
are homogeneous. This raises the need to restrict all integrations to a fiducial cell V, where all 
edges of the cell align with the coordinate axes Xi and with a chosen metric qab- For Bianchi I, 
this metric can be chosen to be flat, and we designate the sides of the cell along the Xi coordinate 
to be Li] thus, the volume of the cell V is given by Vq = L1L2L3. Then, we can write the triads 
and connection components as 



where q is the determinant of the fiducial metric Qab- Throughout this paper, there is assumed 
to be no summation if all repeated indices are covariant or contravariant; the usual Einstein 
summation convention will hold when there are mixed indices. Thus, there is no summation in 
either equation listed in (2.1). The phase space of the reduced symmetry system is coordinatized 
by six variables pi, c*, satisfying the Poisson bracket relation 



Here 7 is the Barbero-Immirzi parameter, a quantization ambiguity which has been calculated 
to be 7 ~ 0.2735 by matching with the Bekenstein-Hawking black hole entropy formula [55]. 
However, there is a freedom to re-scale^ the coordinates independently in Bianchi I as 
under which — a*a;* and — )• (a*) ^e" (no summation on i. Under this re- 
scaling, the variables pi, cP will similarly change - an undesirable feature, since the quantization 

^In most works on LQC, these fiducial inputs are designed as, e.g. °e", to differentiate them from the physical 
analogues el, which includes measurable functions such as the scale factors aiij). However, we do not need this 
distinction, so we simplify our notation accordingly. For more discussion, see [5, 9]. 

^Although we do not address it explicitly here, there is also the issue of the behavior of the system under 
parity transformations [6, 10]. The appropriate choice (when there are no fermions) allows an effective reduction 
in the parameter space m over which the wave function is defined, e.g. only > for all three parameters m,. 





a) 




{c\pj] = ^TTG^V^-H]. 
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of the model should not depend on the fiducial structure chosen. Thus, we normalize the va- 
riables Pi , such that 

Pl = L2L3P1, P2 = L1L3P2, P3 = L1L2P3, 

and 

= Lic\ = L2C^, = L^c?. 

Another way to say this is that the presence of the fiducial volume Vq gives a one-parameter 
ambiguity in the symplectic structure (2.2), which is now fixed by the choice of c*, pj. Therefore, 
the Poisson brackets become 

{c\pj} = S7,G^5]. (2.3) 

If the structure constants C*^ are non-zero - as in the other Bianchi cosmologies - there is the 
need to normalize them accordingly [10]. 

In the full theory, the variables are the fluxes of the triads £"? and the holonomies he , defined 
by a connection A\ along an edge e, where the latter is given by he{A) = V exp{f^Al^Tie°'dt). 
Here, Tj are the generators of the SU(2) Lie algebra, such that TjTj = (l/2)e^jTfc — (l/4)5ijl (these 
are related to the Pauli matrices cjfc as = — ic7fc/2), and I is the unit 2x2 matrix. The fact 
there are no quantum operators corresponding to the connection components c* in LQG is an 
indicator of the diffeomorphism invariance of this theory. To match LQG as much as possible, 
we exponentiate the connection components appropriately to get the holonomy h^j^^ along the 
coordinate axis Xk with an edge of length 5Lk, given by 

4'^ = e^--=cos(^)l + 2sm(^)r.. 

Here (5 is a constant value and the holonomy is given in terms of the almost periodic functions 
exp{iSck) - "almost" because 6 can be any real number, not just an integer. 

We consider a Hilbert space basis of eigenstates |mi, 7712,7713) = \rn) of the momentum ope- 
rator Pi. Here, the values rrii G M may either be the momentum eigenvalues themselves - 
i.e. Pi{rh) = rrii, and therefore have the appropriate physical dimensions - or dimensionless 
quantities, so that the dimensions are carried within the function pi{rn), depending on the 
choice made in a particular physical model. These eigenstates satisfy 

{mi,m2,m3\m[,m'2,m'3) = (5m,,m'/m2,m^,^m3,m?j- 

Because there is a Kronecker delta on the right-hand side, rather than a Dirac delta distribution, 
the wave function is made of a countable number of these orthonormal basis states, in the 
following form 

mi,m2,m3 

mi,r?i2,m3 

where the usual finiteness condition 

I '/^mi,m2,m,3 ' 

mi ,7712, ma 



^ ^ |V'mi,m2,m3| < OO, 



is assumed. This quantization is in-equivalent to the standard Schrodinger form, and thus to 
the usual quantum cosmology using the Wheeler-DeWitt equation as a starting point. When 
the triad and holonomy are promoted to quantum operators, they act on the states as 



Pi\m) = pi{m)\m) , exp{i6c^)\rh) = |mi + fed, 7712, ma) 
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for a constant k and similarly for the other holonomy operators. The function Pi{rn) and k must 
combine to satisfy the Poisson bracket relation (2.3); for Bianchi I, the triad eigenvalue pi{fh) = 
nii and k = —STrGjh. Thus, we see that when the Hamiltonian constraint C = Cgrav + C'matter 
acts on the wave function, since the operators either give functions of fh or shift the basis state 
|m) — )• |m + k6), the equation Cip = gives a difference equation for the coefficients ipmi,m2,m3- 
Up to now, we have considered (5 to be a constant; however, when we consider quantization 
methods in the following section, we will see it is physically more viable to use a functional form 
6 = S{m). 

2.2 Quantization methods 

We consider here only the gravitational piece Cgrav of the Hamiltonian constraint, coming from 
integrating the full constraint over the fiducial cell V, 

C'grav = / NT-Lg-riwdV. 

Jv 

The Hamiltonian ^grav depends on the curvature tensor F^^, so to quantize this tensor, we write 
it in terms of our classical phase space variables, and then promote them to quantum operators. 
In the classical theory, this is done by considering a plaquette - a closed loop whose sides 
are parallel to the coordinate axes and - and defining 

F\ = 2 hm Tr ( ^"'^ ~ ) ojtooi. 



Here, Arn is the area of the plaquette and the holonomy /in-^ is taken as 

where SiLi is the length of the ith edge of the plaquette, measured in the fiducial metric Qab- 
In this scheme, the area of the plaquette is shrunk to zero, so the particular choice of loop is 
irrelevant. 

However, this procedure cannot be done in LQG, due to the discreteness of the spectrum 
for the area operator Ar; instead, the plaquette can only be shrunk until its size reaches the 
minimum area eigenvalue A^p of the area operator, where is the Planck length'^. There 
appears in the literature three separate methods for doing this, with various physical rationales; 
we discuss these here, and in particular their impact on the eventual difference equation. 

• Ho scheme: In this case, the (matrix elements of the) holonomy hj^^'^ is thought of as an 

eigenstate of the area operator Ar j = \pj \ ; the accompanying eigenvalue is set equal to the 
area gap A.ip. In other words, with 

for some constant C, we let 

This gives a constant value for /Uq; this scheme was used by Ashtekar, Pawlowski and 
Singh for the k = isotropic case [5], and Chiou for Bianchi I [34]. The advantage of 



^From the full theory, we have that the dimensionless number A = 47r7v^ although the value A — 2-k'y\/3 is 
sometimes used in earlier works; see the discussion in Appendix B of [1] for further detail. 
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this method is that the difference equation has constant steps, i.e. defined on a lattice 
with even spacing, and therefore is relatively easy to find solutions. However, this scheme 
was realized to be problematic in terms of physical predictions [1]. For instance, as we 
will explain in Section 4, although the singularity is resolved - i.e. the cosmological model 
exhibits a "bounce" if a scalar field is chosen to be the time parameter - the critical 
energy density pcrit at which this occurs is dependent on the momentum of the scalar 
field, with pcrit ~ P^^ ■ Thus, a semi-classical state for the field, with large p^, would 
have a correspondingly low value for pcrit; indeed, there is no lower bound for this density, 
so even densities considered far in the classical regime could be critical densities for the 
appropriate system. 

To mitigate these issues, the constants 5i associated with each edge of the plaquette can be 
generalized into functions 5i{m) of the parameters. Two possible ways to do this are considered 
below; we use here the notation of Chiou [35]. Both of these are dependent on the LQC spin 
network state \p) - parameterized by the triad eigenvalues Pi{m) - and its relation to the 
plaquette chosen. 

• scheme: This particular scheme was first suggested by Chiou [34], and it is inspired 
by the fact that physical areas in LQG are associated with edges between vertices in the 
spin network. Thus, we imagine the side of the fiducial cell V normal to the X3 axis to be 
pierced by N3 edges, each of which is associated with an area and giving a plaquette O12. 
Since for a state \pi,P2,P3) on the kinematical Hilbert space, the area of the face of the 
fiducial cell, orthogonal to the direction Xi, is given by pi, we have [9] 



Since we are assigning each of these edges in the X3 direction to an area parallel to the 
1-2 plane, this gives a plaquette with sides of coordinate size 6^{p)Li and 63{p)L2, as in 
Fig. 1 [35]. The total area of the 1-2 side must be L1L2 in the fiducial metric, so that 



Solving for S3{p), and making the same argument for the 1-3 and 2-3 sides, we have that 
the functions Si{p) in the Bianchi I model are 



In Bojowald, Cartin and Khanna [21], this is the case where the number of vertices is 
proportional to the transverse area. This method is similar to that of the /io scheme, 
in that the difference equation for the Bianchi I model is defined on a constant spacing 
lattice, after the appropriate redefinition of variables [34]. However, like the /Uq scheme, 
working through the effective dynamics of the system shows that its physical predictions 
are dependent on the choice of fiducial cell V - in particular, the densities at which a bounce 
occurs in one of the three spatial directions [35] . 

• ft' scheme: This method was briefly considered by Chiou [34], and developed in detail by 
Ashtekar and Wilson- Ewing [9]. The first step in this scheme is similar to that of the 
method, namely considering the edges passing through a particular face of the fiducial 
cell V. Thus, we again have the relation (2.4) for the relation between the number of 
edges passing through the 1-2 face and the corresponding value of the triad parameter p^. 
However, we now divide this face into plaquettes with lengths 5i{p)Li and 52{p)L2 in the 
fiducial metric as shown in Fig. 1, with the total area equaling L1L2, so that 



(2.4) 



Ns5l{p)LiL2 = L1L2. 




N^6i{p)Li52{p)L2 = L1L2, 
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Figure 1. Illustration of the fl (left picture) and p.' (right) quantization schemes, described in the text. 
The shaded regions in each diagram indicate the area set equal to the minimum area eigenvalue A^p, 
with the corresponding sides labeled by the edge lengths. Figure adapted from Chiou [35]. 



or 



A/2 

If we iterate this for the other two planes, we find for Bianchi I that 

6,{p)=epJ^, 62{p) = £pJ^, 5-M=^P\—- (2.5) 

This is where the number of vertices in a particular direction is proportional to the exten- 
sion in that direction [21]. It can be shown that the physical predictions of this version do 
not depend on the original fiducial cell^, so the critical density is invariant under changes 
in any scalar field used as a time parameter. However, it also results in a difference equa- 
tion for the Hamiltonian constraint where the step sizes are parameter-dependent, so the 
lattice is not made of constant sized steps. From the analogous situation in computational 
physics, we denote this as a "lattice refinement" model. 

Note that the Jl and Ji' schemes are equivalent when there is only a single triad parameter, as 
in the isotropic cases. In the literature, the method is the "original" quantization scheme, 
while Jl' is called the "improved" quantization. Nevertheless, other schemes are possible in prin- 
ciple; see for instance [62], where isotropic embeddings of consistent quantizations of anisotropic 
models are considered. 



3 Suitability of the difference equation 

When one finds the numerical solution of a differential equation, representing the time evolution 
of the system, there is a good deal of freedom in choosing the particular discretization scheme 
for this problem, since small differences in this scheme will be of higher order, and thus drop 
out. Thus, such numerical algorithms are chosen to incorporate a variety of nice features - not 
only matching the original equation in the continuum limit, but also stability under numerical 

''strictly speaking, this lack of dependence on the volume of the fiducial cell V is true only for the effective 
semi-classical equations, but does not hold for the quantum theory. Indeed, looking at holonomy corrections at 
all orders, one finds the expectation values (V') and iV'^) depend on the fiducial cell due to the different scaling 
properties of various pieces of these functions (see Section IV of Chiou and Li [38]). However, this is an expected 
occurrence, as argued by Chiou [35], analogous to the conformal anomaly in quantum field theory; physical input 
on the "correct" size of the fiducial cell may obviate the worry about scaling dependence in the JX scheme. 
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computation. The situation is the opposite for LQC, where the difference equation is fundamen- 
tal, and cannot be altered. However, there is still some ambiguity in the model, dependent on 
the quantization method used to arrive at the final equation; the suitability of these equations 
depends not only on mathematical arguments [40] (such as existence of an inner product in the 
Hilbert space), but also numerical factors, much as the case with differential equations. 

In this section, we consider several issues important in the process of numerically solving the 
Hamiltonian constraint for LQC. First, we must ensure that the difference equation is stable, 
meaning that solutions "close" in the appropriate sense do not exponentially separate as the 
system is iterated by increasing one of the model parameters. This leads to the notion of von 
Neumann stability analysis. Next, the quantum constraint equation must have an appropriate 
semi-classical limit. One must be careful that this limit is unique, which has not been the case 
for all models - the Schwarzschild interior will serve as our example here. Finally, it is helpful 
to consider changes in variables which decrease the difficulty in solving the difference equation. 
A particular problem, coming from the improved quantizations detailed earlier, is that the lattice 
spacings are no longer constant, but depend on the model parameters themselves. In some cases, 
this can be alleviated by choosing a new parameter as a function of the old, where the lattice 
spacing becomes constant up to higher order terms. 

3.1 von Neumann stability analysis 

In the latter half of the 20th century, techniques were developed to numerically compute solutions 
for important differential equations by discretizing the equation and evolving in time on a lattice 
where each point represented a specific space-time location. It was quickly realized that such 
solutions could grow without bound if a bad choice of discretization or time step was made. 
Thus was born the idea of stability analysis, to check the algorithm used and ensure solutions 
remain of finite size. Since the difference equation in LQC models is fixed for a particular 
quantization scheme, the role of stability analysis here is to provide input on the viability of 
the quantization chosen; the assumption is that an unphysical quantum Hamiltonian constraint 
would have several problems indicating its lack of consistency - not only on the basis of its 
predictions, but also from the characteristics of its numerical solutions. This notion has been 
born out for several models displaying many undesirable features, such as the isotropic subspace 
of Bianchi IX [22], the k = \ Friedmann-Robertson- Walker (FRW) model [44], and both locally 
rotationally symmetric (LRS) Bianchi I and the Schwarzschild interior [64]. 

We summarize here the discussion of Bojowald and Date [22]. In their work, they ensure 
the notion of a stable difference equation by checking whether two solutions with "nearby" 
initial data will remain close to each other, rather than diverging. For simplicity, we start with 
a difference equation for a sequence V'm of one parameter m, or 

k 

2^ Ai{m)^„,+i = 0. (3.1) 

i=—k 

The analogous analysis can be done for multiple parameters as well. We look at small regions 
in the parameter space, so that the coefficients Ai{m) of the equation can be approximated by 
constant values Ai. In general, solutions of difference equations with constant coefficients are 
found by writing down the characteristic equation 

k 

^ AiX' = 0, (3.2) 

i=—k 

and solving for the roots Xa [42]. If there are no repeated roots, the solutions are linear combi- 
nations of the form ipm = /3iA]" + • • • + /32fc+i-^2fe+i' with complex constants For a single 



8 



D. Brizuela, D. Cartin and G. Khanna 



repeated root Xa of multiplicity da, the basis is modified to {A^", A2^, . . . , A™, mX^, . . . , m'^'"~^A™, 
• • • '^^21+1}' similarly for more than one repeated root. As shown by Bojowald and Date, 
the physical requirement that the difference equation has the appropriate classical limit means 
that A = 1 is a double root of the characteristic equation. 

The need for a stable difference equation means that two sequences Vm and ip'^, star- 
ting at a parameter value M with initial data {ipM,ipM+i, ■ ■ ■ , V'Af+2fe-i} and {V'm' V'm+h • • ■ ' 
V'A/+2fc-i}' ^^^^ remain close together if the initial data is close as well. These sets of values 
and the difference equation (3.1) allow one to solve for ^M+2k, and all subsequent values for 
m > M + 2k. If we think of each set of initial data as vectors S{M) and S'{M), respec- 
tively, in C^*"' and use the Euclidean norm || • ||, then the initial data are close to each other if 
US' — 5" 1 1 < 6m- Then we require that the sequence 61pm = i^m — 'ip'm remains small in size as well; 
this places conditions on the roots of the characteristic equation (3.2). We see this as follows. 
For small ranges of the parameter 5m, ^ipm satisfies the same constant coefficient difference 
equation (3.1) that the individual sequences ipm and •i/'m do. Thus, we can write 5ipm as a linear 
combination of the basis elements as 

2k da-1 

£1=1 fa=0 

where a ranges over the separate characteristic roots, and Ca,ra are constant coefficients. For 
a given sequence tpm and set of initial data S{M), we chose a nearby initial data set S'{M) such 
that only one of the values Ca,ra is non-zero. Then as the sequence 5ipm goes from parameter 
values M to M + AM, we have the norm of the difference between the two sequences goes as 

<^Af+AM — |Aa|^*^ ^1 + 1 ^ ^M- 

In order to keep 6ipM+AM within the same bound on the norm as 6^pM, we must have all the 
roots |Aa| < 1. 

In Appendix A, we examine the stability of two versions of the k = 1 FRW model, that of 
Green and Unruh [44], the other by Ashtekar, Pawlowski, Singh and Vandersloot [7]. Although 
the Green and Unruh model is unsatisfactory due to its unphysical predictions, the stability 
analysis shows that the problem lies with the quantization method and the resulting Hamiltonian 
constraint, an issue repaired in the work of Ashtekar et al. 



3.2 Semi-classical limits of the difference equation 

When using these discrete equations to obtain predictions for semi-classical effects, one must 
be careful about how the limit is taken. We assume that the model difference equation is 
essentially continuous in the semi-classical limit, in order to match with the existing theory of 
general relativity. In other words, we should obtain the WdW equation as the limiting case of 
the difference equation. For LQC, the (dimensionless) quantum parameter m is related to the 
continuous variable p{m) by 

p(m) = jim, 

where 7 is again the Immirzi parameter, and £ is proportional to the square of the Planck 
length. Since p{m) is the eigenvalue of the triad operator, it is related to the scale factor func- 
tions appearing in the classical metric. One possible way to take the semi-classical limit is letting 
7 — )• 0; however, since 7 is a small but fixed constant, we must assume that 7 — )• 0, m — )• cxd such 
that p{m) remains constant. This is what is called the "pre-classical limit" in Bojowald and 
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Date [22]. This respects the physical input that 7 represents. Then, we can approximate the dis- 
crete evolution equation (3.1) with a differential equation for a smooth wave function %l){p{m)). 
This can be obtained by writing the wave function V'm and the coefficients ^4^(7) in a power 
series in 7. Putting these expansions back into the original discrete relation, we find a con- 
tinuous differential equation as the leading order contribution. For LQC, the WdW equation 
should emerge as the 0(7^^) term, with the higher order 7 terms corresponding to quantum 
corrections, and all terms with inverse powers of 7 should cancel out. The conditions to obtain 
such a physically reasonable situation are derived by Bojowald and Date [22]. 

In some sense, taking the limit as described above is an example of finding the semi-classical 
equation first, then solving the resulting WdW equation for the wave function. However, we 
could solve for the wave function first without appealing to the physical nature of the parame- 
ters - using generating function techniques [30, 31, 32, 33] or other methods - and then pick out 
those solutions that have the appropriate behavior. This corresponds to taking the asymptotic 
limit m — )• 00 without reference to the parameter 7. This is the "pre-classical limit" as it has 
been used by Cartin and Khanna [31, 32, 33]. The leading order term in a m power series 
will remain a difference equation. Generally, solutions of difference equations can have many 
undesirable properties, such as oscillation in sign with each succeeding step (i.e. ■0m — (—A)™' 
for a positive constant A in the large m limit) or unstable solutions that increase without bound. 
Some control over these solutions can be obtained by choosing the proper initial data. However, 
it is possible that there is not enough free parameters to avoid unphysical behavior. 

As an example, we look at the original quantization of the Schwarzschild interior [2, 58]. If 
we write it in terms of the two triad components ph, Pc (because of spherical symmetry), the 
Hamiltonian constraint is of the form 



+ 



-ffpS(^^Pc + 2jilS + [ij{pb + ^il5,pc + 2jfp6) - ^'(Pfe - 7fpS,Pc + 27^5)] 

^Jpc + lilS - yjpc- 74^^) [iPb + l^l^MPb + 2^fp6,pc) - 2(1 + 2j^5^)pbilj{pb,Pc) 
+ {pb-lilS)i^iPb-2jfp6,p,)] 
+ ^el6(^^Jpc-2j£l6 + [i;{pb - ^fp5,p, - 27^^) - V(Pfe + lilS,pc + 27^5)] = 0. 

For most of the terms and coefficients of the difference equation, taking the Pb,Pc 00 limit is 
equivalent to that of 7^p5 — )• 0. However, there is one problematic term in the equation, namely, 

T{pb,Pc) = -2Upc + 7£l6 - ^pc - lip) (1 + 2jH^)pbij{pb,Pc). 



Specifically, if we look at the quantum gravity limit j6 — )■ 0, but allowing pb, Pc to remain the 
same size by requiring mi,m2 — )• 00 accordingly, then 

lim T{pb,pc) = -{2j6ij 



while in the large momentum limit, then we have 



2x2\ Pb 



lim Tipb,Pc) = -(27^4) (1 + 27''5- 



Pt,Pc-^O0 ^p^ 



The difference in the limits means this model has the undesirable feature that the nature of the 
semi-classical limit changes, depending on how this limit is taken. Thus, the best case is to find 
a Hamiltonian constraint and difference equation where the semi- classical regime is the same, 
regardless of which limit is used to reach it. 
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3.3 Lattice refinement and model reparametrization 

The advent of these improved quantization schemes unfortunately leads to difficulties in sol- 
ving the difference equation, which now involves lattice refinement. As we saw above, the 
and p,' schemes both involve lattice steps that are functions of the model parameters m. If 
we have a difference equation of two parameters, say mi and m2, then generic terms in the 
Hamiltonian constraint for the wave function 'ip{mi,m2) will be of the form f{mi,m2)ip[mi it 
(5i (mi, m2), m2 ±52(^^1 ) "12)] where the step functions 6i{mi,m2) are non-linear functions of the 
parameters. These difference equations are more difficult to solve, since unlike relations with 
fixed steps, the equation for a given choice {m} depends on values of ip(rn) not calculated from 
any other equations. An example of this is the difference equation arising from a lattice-refined 
Schwarzschild model [21], where 



di(mi,m2)= . 02{mi,m2) = . (3.3) 

Vimal mi 

Note that 5o is a fixed parameter related to the minimum area eigenvalue (i.e. 60 ~ VA)- In 
this case, the Hamiltonian constraint is of the form 

C+(mi) [tp{mi + 260/ \/|m2|, m2 + 25o\/k"2|/"T-i) - tp{mi - 25q/ \/\m2\, 

7112 + 25Q\/\m2\/mi)\ + C°(mj) [(mi + 25q/ \/\rn^\)'ip{'mi + A5q/ \/\rn^\, m2) 
- 2(1 + 27^(^o/"^2)miV'(mi, m2) + (mi - 26q/ \/\rm\)'4){mi - A5o/ \/|m^, m2)] 
+ C^(mj)[-0(mi - 25o/ ^/\m^\,m2 - 25Q\/\m^\/mi) - V'(mi + 25q/ \/\rri^\, 
m2 - 2(^0 Vl "12 1 /mi)] = 0. 

To solve these types of lattice refined equations, we consider a change in the parameters used 
in the model, where the equation in the new variables now has constant size steps [21, 59]. This 
is possible always for one-parameter models, i.e. isotropic cosmologies, but may not be feasible 
for equations with multiple parameters. This can be seen as follows. If the wave function 
depends on only a single variable m, then with non-constant steps 5(m), we have terms in the 
difference equation of the form ip[m ^ 5{m)\. We assume the step function 5{m) is linear in a 
parameter 5o which we use as our lattice spacing. To get constant step sizes, we must write the 
wave function tp{m) in terms of a new function N(m), such that 

N[m±6{m)] = N{m)±6o + 0{6^). (3.4) 

This would imply that 

ip[N{'m ± 6{m))] ~ iplNim) ± 5o], 

so the values N{m) form a lattice with a constant spacing 5q. Taking the Taylor series expansion 
of the left-hand side of (3.4) gives that this condition requires 

dN{m) 

— d{m) = 60. (3.5) 

am 

If the function N{m) is such that higher order terms in (3.4) correspond to higher order in h as 
well, then they can be dispensed with using the appropriate factor ordering of the Hamiltonian 
constraint operator [21]. Nelson and Sakellariadou showed the issue of factor ordering is inti- 
mately related to that of lattice refinement by showing that, for the flat isotropic model, there 
is a unique choice of quantization scheme (the "improved" version) where the factor ordering 
ambiguities disappear at the semi-classical level [60]. 
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An example of this reparametrization procedure is seen for the flat isotropic model [6] . The 
authors of that work use a geometric argument - considering integral curves of operators related 
to the holonomy operator - to find the new function N{m) such that ip{N(m)) has a constant 
size lattice. Here we show the same result is obtained by the method described above. In 
particular, in this isotropic model, the need to obtain valid physical predictions (much like for 
the Bianchi 1 model) makes it natural to have a step function 5{m) = 2/{2>K^/rn)^ where K 
is a purely numerical constant. Solving the equation (3.5) for N{m) with 6q = 1 gives N{m) 
is proportional to the volume eigenvalue associated with m, that is, v{m) = A'sgn(m)|m|^/^. 
Therefore, if we redefine the wave function as tlj{v) = ip{m), then the action of the operator 
with the improved quantization function p,{m) gives equidistance steps. 

This method becomes problematic for multiple parameter models. As we have seen with the 
functions 6i{m) - for Bianchi I, given for the jl' scheme in (2.5), while in the Schwarzschild 
interior, shown in (3.3) - the step functions Si can depend on all parameters mj. Suppose for 
specificity we have two parameters. Then, we seek to find two independent functions N (mi, 1712) 
such that 

N[mi ± ki6i{mi,m2),m2 ± A;2(52(mi, 7712)] = N{mi,m2) + k{ki, k2)6o + 0(^5q), 

for some constants ki, and k a function of fcj. This requires that 

dN dN 

(5i(mi,m2) + ^ d2{mi,m2) = kSo. 

omi 01712 

For this is to be true for all ±ki, ±k2, we need specifically for the Schwarzschild case that 
2 dN ^ 2./\^\ dN ^ 



y^\rn^\dmi ' nii dm2 

for constants Q; these equations only have one independent solution, namely N{mi,m2) ~ 
miy^|m2|, i.e. the eigenvalue of the volume operator for the model. Therefore we see that it is 
possible to have constant step sizes in one of the two parameters, but not in both [21]. This 
simplification, however, is useful in putting the difference equation into a form more tractable 
for numerical computation [46]^. 

Beyond this analytical technique to place the original difference equation into a more tractable 
form, there are also various methods of numerically solving the equation along the lines of the 
Taylor expansions outlined above. One possibility is to interpolate the solution between lattice 
points already calculated, using a least-squares fit [65] (more generally, the wave function is 
expanded in a Taylor series of the parameters around these previously computed values [59]) 
and then use this interpolant to compute the values needed to step the evolution forward. 



4 Examples of LQC models 

In this section we will mainly restrict our attention to those studies that have faced the purely 
quantum difference equation via numerical techniques. This has been performed for a number 

^Solving a partial difference equation with non-uniform stepping is not straightforward, i.e., a naive iterative 
approach towards computing the solution simply does not work. The reason is that, due to the non-uniform 
stepping, computing the future value of the solution requires one to know past of values of the same, that happen 
to be different from the ones computed in previous iterations! However, one can overcome this issue through the 
use of a local interpolation formula for the solution, that enables one to compute precisely those values of the 
solution that are actually needed, thus making it possible to propagate the evolution forward [65]. Now, if one is 
able to transform the equation into a form in which one of the variables has a constant step-size, then the local 
interpolant is not required for that variable, and one can proceed with the computation through a simple stepping 
procedure on that variable. This simplifies the computation considerably [46]. 
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of homogeneous models: FRW with spatially flat k = 0, closed k = 1 and open k = — 1 
topologies, with and without the presence of a cosmological constant and, less systematically, 
for anisotropic Bianchi I models and Schwarzschild interior. The quantization of Bianchi II [10] 
and Bianchi IX [75] models have also satisfactorily been carried out in LQC, but the numerical 
analysis is still missing. In a last subsection we will mention several studies where the effective 
analysis has been done, since we consider this is a very important research area of LQC. On 
the one hand, many of the main features of the LQC evolution has been first derived via this 
effective treatment. On the other hand, this effective analysis have been applied to a wide variety 
of scenarios, such as inflation as well as to some anisotropic and inhomogeneous model, which 
permitted to extract physical information through systematic numerical evolutions that would be 
difficult to perform in the full quantum setting. Of course, since to a certain extent, heuristics is 
used within these approaches, care is needed when interpreting the validity of the physical results. 

4.1 Flat Priedmann— Robertson— Walker 

Because of its simplicity and the fact that it describes to a very good approximation the evolution 
of our universe, the FRW solution has been the most studied model in LQC literature. In fact, 
the first positive results regarding the singularity resolution in the context of LQC were due 
to Bojowald [14, 16] for this model. The singularity is indeed resolved not simply because the 
evolution equation is discrete and, hence, the state jumps over the zero point. In fact, the 
evolution through the classical singular point is also well defined. Therefore, any initial data 
can be in principle deterministically evolved through the singularity. But, in order to know 
precisely what is "on the other side of the singularity" and what is the precise behavior of the 
physical system, one needs to construct Dirac observables and solve the difference evolution 
equation. Here is where one has to resort to numerical methods. The seminal work in this 
aspect was performed by Ashtekar, Pawlowski, and Singh (APS) [5]. Let us comment more 
precisely about their approach. In their work, the spatially flat FRW model was analyzed with 
a matter source of massless scalar field 0. This model is particularly interesting for the study of 
singularity resolution since classically a (big-bang) singularity is unavoidable. The scalar field 
is a monotonic function of the cosmological time, so it can itself be interpreted as an internal 
clock instead of the scale factor a. On the one hand, this is a widely used approach in cosmology 
to describe the evolution of the system without making reference to any fiducial structure or 
coordinates. On the other hand, in the closed (k = 1) model, although the universe undergoes 
a recollapse, the scalar field (j), contrary to the scale factor a, is still a monotonic function and 
can be interpreted as a time variable. 

In this model, the phase space is coordinatized by the gravitational degree of freedom {c,p), 
with Poisson bracket {c,p} = SirG/S, by the scalar field (j) ^-nd by its conjugate momentum 
with bracket {(j),P(j>} = 1. The only constraint left is the Hamiltonian that, in terms of these 
variables, takes the the form, 

C = N[Cgra.v + C'matter)) 




For this system, a complete set of classical Dirac observables is provided by the constant of 
motion and by p|<^o- "^^^ latter determines the value of p at the instant of time at which 
(j) = (pQ. This set is complete because they univocally define a trajectory on the phase space. 
Therefore, in order to extract physical information of the quantum system, the evolution of these 
observables must be analyzed. 

The quantization of the geometric degrees of freedom is performed by following the methods 
of full LQG and therefore, as has been explained in Section 2.1, making use of holonomies A'^^ 




with 




and 




.|3/2- 
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and triads p as fundamental variables. In this way, the gravitational kinematical Hilbert space 
is 'H^^ = L'^{M.B,dfiB), being the Bohr compactification of the real line with the Bohr 
measure d/j^B- This Hilbert space is properly spanned by the eigenfunctions |^) of the triad 
operator: p\fi) = Sirjip^fi/dlii). On the other hand, the matter degree of freedom (j) is quantized 
via a standard Schrodinger representation on a Hilbert space 'HJJJ^** := L'^{'M., dcj)). Here, the 
operator associated to the conjugate momentum of the matter field p,^ acts as a derivative 
operator p^p = —iKd^j,. 

In this first study by Ashtekar, Pawlowski and Singh, the lapse was taken equal to unity 
and the //q scheme was used. The physical states were then defined as those that are 

annihilated by the Hamiltonian constraint. On the full kinematical Hilbert space, given by the 
product between the geometric and matter spaces L^(Mb, B{p)dnB) ® d(j)), the constraint 

can be written as an evolution equation in terms of the internal emergent time 0, 

dp! = -eo^-, (4.1) 

where the difference operator is given by, 

eo^(/i, 0) := [B{p.)]-'{f{p)^{^^ + 4;Uo, <P) - [f{p.) + /(^ - 4^o)]^(/i, <A) 

+ f{pi-A^l^)^{^l-^^lQA)], (4.2) 

for certain functions f{p) and B{p). In particular, B{p) is proportional to the eigenvalues of 

the operator l/|pp/2, which is the only geometric factor that appears on the matter part of the 
Hamiltonian. Due to the difference operator, superselection sectors appear. Let us denote by 
the lattice of points {e + 4/io} on the n axis. Therefore, the subspaces Tie G T-lkm, defined as 
those states with support only on Cg, are left invariant under the action of the operator ©o and 
the Dirac observables. 

The positive-frequency solutions to this equation form the physical Hilbert space and they 
can be written in an integral form as, 

/oo 
dk^ik)e)!^e'''^''^'f'. (4.3) 
-oo 

In this expression, ^(k) is an arbitrary function, uj{k)'^ = 7rG(16fc^ + l)/3, whereas e^*^ is 
a symmetric linear combination of eigenfunctions of the operator Oq with eigenvalue k. These 
combinations have support on the lattice and are, by construction, invariant under the triad 
orientation reversal operation. The physical inner product for a given value of (p = (po is, 

(^'l|^'2)e= Yl ^(/«)^l(/">o)^2(Ai,(/'o), 

and the action of the Dirac observables is given as. 



|//|^3*(^,</)) = e*^('^-^oVl^(/^,'/'o), and p^^(/x, </.) = -i/i- 

In this way, there are two possible methods in order to solve the constraint equation. On the 
one hand, the constraint can be understood as an evolution equation in 4>. Hence, choosing initial 
data "^{fi, 4>o) at a given value of the scalar field (j) = (po, the physical state is obtained by evolving 
them via equation (4.1). On the other hand, one could obtain the eigenfunctions of the opera- 
tor 00, construct the appropriate symmetric linear combinations e^*\ pick up a function ^(fc), 
and perform the integral on the right-hand side of equation (4.3). The advantage of the first 
method is that it is not necessary to deal with the eigenfunctions of Oq but, from a numerical 
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point of view, it is more difficult since a large number of differential equations must be solved. 
The choice of the free data in each procedure (either as initial data or as the free function ^) will 
introduce the physical input on the system. Obviously, the most interesting set of initial data 
is given by those that describe a large classical universe at late times, since this is the current 
state we observe in our universe. In [5] these both procedures were followed to numerically solve 
the behavior of the system. Let us comment on those with some more detail. 

For the evolution procedure, three different methods were used to pick up initial conditions. 
In the first one a Gaussian (with respect to the Wheeler-DeWitt inner product) peaked on 
a classical universe was chosen. In the other two, an analytical solution of the Wheeler-DeWitt 
equation, that is semi-classical at late times, was used. Since equation (4.1) represents an 
infinite number of coupled differential equation, one needs to restrict to a finite domain in order 
to perform the numerical evolution. The boundary was chosen far away from the peak of the 
initial wave packet, and purely outgoing boundary data were chosen there. In addition, for the 
discretization of the time derivatives of (4.1), the fourth-order adaptive Runge-Kutta method 
was applied. Regarding the direct evaluation of the integral (4.3), the function ^(k) was chosen 
as a Gaussian (with the usual measure) peaked on certain value k* (<C — 1) so that it describes 
a semi-classical state at a time (pQ. The eigenfunctions of ©o were numerically obtained and the 
integral (4.3) was evaluated using fast Fourier transform. 

Even if, obviously, the considered different initial data sets led to different quantitative 
results, the generic behavior and the main qualitative results were robust independently of the 
method of choosing the initial data or the value of the free parameters within those methods. 
Evolving backwards, the expectation value of the observables follow the trajectory of their 
classical counterparts until a certain minimum value. At this point they depart completely 
from the classical behavior and they bounce, keeping always a finite value, and connecting an 
expanding with a collapsing branch of the universe. In addition, the state is sharply peaked 
during the whole evolution. In this way, the classical singularity is resolved and the LQC 
evolution equation connects two classical branches of the universe through a quantum phase. 
The bounce happens because the system reaches a point where quantum gravity effects are 
appreciable and make the gravity repulsive. In principle, one would assume that this must 
happen for high values of the matter density, since we do not observe any quantum gravity in 
our daily experience. Indeed, this fact is one of the drawbacks of the analysis presented in [5] 
since the matter density at the bounce happens to be inversely proportional to the expectation 
value {Pff,), as discussed in Section 2.2. Therefore, one could pick up a value of {prj,) so that 
the bounce would occur at very low values of the matter density (and, hence, curvatures). This 
drawback was fixed in [6] by introducing the so-called improved dynamics. 

In this analysis, the fl quantization scheme was used. Note that in this homogeneous scenario 
the p, and p,' schemes presented in Section 2.2 coincide. The mathematical structure of the 
formalism is essentially the same as for the case explained above, but physically the main flaw 
of that approach was cured. Within the improved-dynamics scheme, for all the considered cases, 
it was obtained that in a backward evolution the expectation values of the Dirac observables 
remain peaked on their classical trajectories (corresponding to an expanding branch of the 
universe) until the matter density reaches a critical value of of /Jcrit = 0. 82 ppi, where the system 
undergoes a bounce. After that, the energy density starts decreasing and, once its ratio with /Jcrit 
is small, the system evolves again peaked on a (contracting) classical trajectory. As already 
mentioned, the value pcrit is the same for all analyzed simulations. In fact, from the effective 
analysis of this system, it is possible to write down a modified Friedmann equation which 
provides an analytical expression for this object: pcrit = 3/ (Wir'^j^G'^h). This value is in very 
good agreement with the one found in the full quantum simulations. 

Regarding the mathematical issues of these second analysis [6] , in order to label the eigenfunc- 
tions of the triad operator, the real parameter is replaced by the affine one v := Jsgn(|u)|/x|^/^, 
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with J = 2v2/ (3\/ 3v3). In this representation, the action of the components of the holonomies, 
g«^/2^ simply shift in a constant step the states, 

e«A/2|i;) = \v + 1). 

This basis is best suited to the volume operator (associated to the elementary cell): 

The quantized Hamiltonian constraint, written as an evolution equation, takes the same form 
as in equation (4.1), with the same structure for the operator, 

Go^(t;, (A) := [C{v)]-^{g{v)^{v + 4) - [g{v) + g{v - 4)]^{v, 0) + g{v - A)^{v - 4, 0)}, 

but with functions g{v) and C{v), differing from B{n) and /(/i) that correspond to Qq (4.2). 
The general (positive-frequency) solution to this equation can be written in a form completely 
similar to (4.3), where now one should consider eigenfunctions of the new difference operator Qq. 
Finally, the numerical techniques to perform both the evaluation of that integral, as well as the 
direct resolution of the evolution equation, are completely analogous to the previous case. 

Here a brief comment about the work by Laguna [49] is in order, which has been a source of 
some confusion (for further discussion, see Section VIB of [7] ) . The paper by Laguna first deals 
with the evolution of wave packets under the LQC Hamiltonian constraint for this isotropic 
model, and then compares it to the WdW equation for that system. Finding a dispersion 
equation for plane wave solutions of these models, it is shown that the group velocity goes 
to zero in the LQC model - and thus the quantum bounce occurs as the contraction of the 
spatial volume becomes an expansion - precisely at the critical maximum energy density pcrit 
shown previously. However, Laguna then goes on to consider modifications of the original LQC 
equation, which has the form of the shallow water equation, that would eliminate the bounce. In 
this context, of interest to those working in the field of computational physics, he essentially adds 
extra terms to the equation, none of which are motivated from the original LQC quantization 
scheme. In this sense, the elimination of the bounce as a spurious effect is not applicable to 
LQC directly, since it involves going outside the quantization method used to arrive at the 
Hamiltonian constraint operator. 

Finally, even if this is the simplest model of LQC, there is still no complete consensus on 
the precise details of the quantization. Even inside the framework of the improved dynamics 
(choosing the p, scheme) there is freedom in the construction of the quantum Hamiltonian since, 
as is usual in any quantization process, in order to represent the classical constraint it is possible 
to choose different factor orderings for the (non-commuting) operators. In addition, the fact 
that one could choose arbitrarily the lapse function also introduces another ambiguity and gives 
rise to different quantization schemes and, hence, to different discrete operators @. In a recent 
paper [56] four of these prescriptions have been analyzed and compared. The first prescription 
is the one due to APS and has already been presented above. In this case, the lapse is taken 
equal to unity and the quantum constraint is symmetrized with respect to holonomy operators. 
In [3] , a simpler version of this scheme was presented in the sense that it allows to carry out the 
complete analytical resolution of the system. This is usually known as solvable Loop Quantum 
Cosmology (sLQC). Its key feature is that the lapse function is chosen to be V/{8'itG), what 
produces a quantum Hamiltonian constraint free of inverse of triads (volume) operator. The 
origin of the third analyzed prescription lies on the analysis of Bianchi I models [9] . As in the 
APS case, the lapse is also taken equal to one but, contrary to the previous prescriptions, the 
operator sgn(z;) appears in the quantum constraint. This prescription is usually denoted as 
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MMO, due to the names of the authors of the article [51], where it was first introduced. Finally, 
in the same paper [56] a new prescription (sMMO) is introduced that combines different aspects 
of the sLQC and MMO ones. 

The system has been analyzed for all the mentioned prescriptions following the route of the 
construction of eigenfunctions of the corresponding discrete operator O. From a numerical point 
of view, it was found that the MMO and sMMO prescriptions are much more efficient. This is 
due to the fact that, for these cases, the spectrum of O is non-degenerate: the state is required 
to be symmetric with respect to a change /i — )■ — and, hence, the mentioned eigenfunctions are 
determined once their value at a given lattice point is fixed. In addition, nontrivial differences 
between prescriptions have been found, although they might be non observable in practice. For 
instance, three observables (the volume, the Hubble parameter, and the energy density of the 
scalar field) have been analyzed and, in all the cases, the difference between their expectation 
values for different schemes turn out to be smaller than their corresponding dispersion. 

4.2 Closed Priedmann— Robertson— Walker 

Quantum gravity effects are supposed to modify the behavior of the system only at Planck 
scales in such a way that the classical singularity is resolved. In addition, they are expected to 
be suppressed at low curvatures in order to recover the classical trajectory. Since classically the 
universe undergoes a recollapse, the closed FRW model represents a more restrictive testbed 
model for LQC methods than the flat case: this recollapse should not be affected by quantum 
geometry modifications. One of the first proposals for the quantization of this model was 
presented in [16] in combination with the flat (k = 0) model. The problem, as has been explained 
above, is that this quantization was proven to be unstable [22] for the closed case. Nevertheless, 
since this isotropic model can be understood as a particular case of Bianchi IX [28] , Green and 
Unruh [44] made use of techniques developed for the Blanchi IX model [18] in order to construct a 
quantum hamiltonian constraint within the /Uq scheme. In this analysis, the numerical resolution 
of this constraint was also obtained. Even if an inner product and observables, that would 
provide a neat physical interpretation of the system, were missing, the divergent behavior of the 
solutions at large scales already indicated that the classical recollapse could not recovered. 

In [7], the fi scheme was used to quantize the closed model. In this article some important 
developments were made that shed light on the concerns explained above. Physically, the most 
important result was that, due to the quantum effects, the classical (big-bang and big-crunch) 
singularities are both resolved and replaced by a quantum bounce. This led to a cyclic universe 
picture, where the quantum wave function evolves through infinitely many classical cycles. The 
value of the matter density at the bounce is in agreement with the critical value pcritj intro- 
duced above for the flat model, (provided that the volume of the universe reaches a macroscopic 
size). More importantly, the state remains sharply peaked on the classical trajectory (within 
the corresponding cycle) and, thus, undergoes the recollapse at a maximum volume that is in 
complete concordance with the predictions of the classical theory. Therefore, the infrared limit 
of the theory turns out to be right. From a technical point of view, the analysis of [7] is formu- 
lated in terms of connections. Contrary to previous quantizations which, for technical reasons, 
constructed holonomies in terms of the extrinsic curvature, in this approach the connection was 
used. This procedure mimics better the approach followed in full LQG. 

The overall structure and analysis is completely analogous to the flat case [6] . The kinematical 
Hilbert spaces are again given by ?^kin- The Hamiltonian constraint can be written as, 

dl^{v,<t>) = -Qi^{v,4>), 

with a discrete operator @i"^{v,(j)) := @o'^{v,(j)) + il{v,fi,ro)'${v,(j)), so that the scalar field cp 
serves as emergent time. The function depends on the radius of the 3-sphere ro with respect 
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to the fiducial metric. The flat k = case is recovered by letting tq tend to infinity, which 
makes the function vanishing. This discrete operator also produces superselection sectors. 
The key difference with the Go operator is that the eigenvalues of 01 are discrete (see [70] for 
an analytical proof). More precisely, eigenfunctions exist for all values of to but, in order to be 
normalizable, they must decay exponentially in both limits v — >• ±00. And this only happens 
for certain discrete values of to. Numerically these eigenvalues were found following a bisection 
method (taking advantage of the fact that the normalizable eigenfunction is always a critical 
solution between functions exponentially diverging to ±00 as v — )• 00). 

4.3 Open Friedmann— Robertson— Walker 

The study of the open k = — 1 isotropic case has been carried out in [74] within the /x scheme. The 
discrete evolution operator 0_i is constructed and the scalar field is understood as emergent 
time. Following the procedure introduced in [5, 6], eigenfunctions of the 0_i operator were 
obtained and then Fourier transformed, in the analogous equations to (4.3), in order to get the 
physical states. The qualitative behavior of the system is completely similar to the plane k = 
and closed k = 1 models described above. Regarding the energy density of the scalar field at 
the bounce, an analytical expression is obtained by effective means, which coincides with the 
results of the full quantum system. This expression depends on the momentum of the scalar 
field but it tends to pcrit for large values of p^. 

In [69] some issues, like the choice of the loop to define the holonomies, were clarified, but 
this analysis is not yet totally satisfactory. First, the defined discrete operator is not essentially 
self-adjoint. Also, as in the first analysis of the closed model, the extrinsic curvature, instead of 
the connection, is used to construct the holonomies. This route is taken in order to avoid the 
technical difficulty of constructing an appropriate loop. In principle this last drawback could be 
cured following the proposal presented in [10] for the Bianchi II model, but a complete analysis 
is still missing. 

4.4 Cosmological constant 

The inclusion of the cosmological constant has only been considered for FRW spatially flat 
models with a matter content of a massless scalar field that, as in previous models, also serves 
as internal time. In particular, in Appendix A of [6], the main techniques for the fl scheme 
were already put forward and some numerical evolutions of this system were presented. In all 
cases the the classical singularity happened to be replaced by a quantum bounce and, contrary 
to previous treatments [11, 63], the correct semi-classical behavior of the system was recovered. 
Moreover, the total energy density (defined as the sum between the energy density corresponding 
to the matter and to the cosmological constant: p := + A/(87rG)) at the bounce was found 
to be the same as in the A = cases and, hence, independent of the value of the cosmological 
constant. 

Already at a classical level, the behavior of the system is completely different depending on 
the sign of the cosmological constant. Classically, the universe with A < begins at a big bang 
singularity at a "time" (j) = —00, expands until the total energy density reaches a minimum 
(in fact it vanishes) where, even if we are in the spatially flat case, undergoes a recollapse and 
finishes at a big-crunch at (p = 00. In [12] a concrete and detail analysis of the quantization 
of this model was carried out and the results were similar to the closed k = 1 FRW with 
a vanishing A. The spectrum of the corresponding discrete operator turns out to be discrete 
and, therefore, its eigenvalues are isolated. Hence, a similar bisection method to the k = 1 case 
was applied. The quantum system also resolves both initial and final singularities, replacing 
them by quantum bounces happening at a value p = 0.82ppi, which results in a cyclic universe 
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model. In addition, the classical recollapse was also recovered in complete agreement with the 
classical results. 

In the classical A > model, after the big-bang, the scale factor of the universe reaches 
infinity at a finite value oi 4> = 4>o. Therefore, contrary to the A < model, in this case the 
values of the scalar field are constrained to a finite interval. In [47] Kaminski and Pawlowski 
constructed the discrete evolution operator making use of the three quantization prescription 
commented above (APS, sLQC, and MMO) and studied its possible self-adjoint extensions. 
The properties of this operator happened to be completely different depending on the value of 
the cosmological constant. For < A < Acrit, where Acrit = SvrGpcrit is the value for which the 
energy density of the cosmological constant equals the maximum total density pcritj the operator 
admits many extensions; whereas for A > Acrit, the operator turns out to be essentially self- 
adjoint. Observations favor a positive but small value of the cosmological constant A <C Acrit and 
the different self-adjoint extensions could give rise to non-equivalent unitary evolutions. Even 
so, a recent numerical analysis [4] shows that for the relevant states (which are initially sharply 
peaked at a classical point of the phase space) the evolution is quite insensitive to the chosen 
extension. Furthermore, the classical initial singularity is resolved. The self-adjoint extensions 
permit to follow unitarily the evolution in (p after the classical divergence at (po and give rise to 
a recollapse. Thus, a cyclic universe picture is again obtained. 

In order to finalize this section, we would like to briefly comment on another approach that 
has been proposed in order to include the cosmological constant in LQC through the so-called 
unimodular gravity [39]. At classical level, unimodular gravity is equivalent to general relativity 
with the only difference that allows for a dynamical cosmological constant. In fact, the momen- 
tum conjugate to the cosmological constant can serve as emergent time and it seems that the 
quantum theory one obtains is simpler than in the cases with the scalar field clock. Remarkably, 
the main results of singularity resolution and the critical value of the energy density coincide 
for both cases. 

4.5 Bianchi I 

Bianchi I model is the simplest non-isotropic model one could consider. Even though, most of 
the analysis for this model, as a test case for anisotropic cosmologies, has been performed at the 
level of the effective dynamics rather than solutions of the Hamiltonian constraint difference 
equation. In part, this is due to the greater complexity of the constraint, with a lattice having 
non-constant steps in the fl' quantization scheme. Complete quantizations have been performed 
both for a massless scalar field matter content [34, 9] and for vacuum [52, 54]. In this section, 
we will just mention these two quantizations and the only numerical analysis [53], even if not 
completely systematic, that has been carried out to analyze the full quantum evolution of this 
model in the fl' scheme. 

As in the previous cases, the scalar field serves as an internal clock and can be treated as 
an emergent time. Making use of this fact, the quantization of this system was first carried 
out in [34] for both /io and fi schemes. Nonetheless, this approach suffered from several draw- 
backs, like the incorrect infrared limit and dependence of the physical results on the choice of 
cell [71], that were fixed by considering the p,' quantization [9]. The particular form of the 
Hamiltonian constraint given in [9] is non-symmetric in its three parameters due to a variable 
transformation chosen by the authors; instead of using all three triad parameters pi to serve 
as variables in the difference equation, one of the three is transformed (as discussed earlier in 
Section 3.3) to a volume v. Thus, the model is given in terms of two directional triad parameters 
and the total volume. In addition, it was shown that the difference equation defined by that 
Hamiltonian is unstable to growing modes [61], although it is unclear whether these modes are 
eliminated by the inner product. An open area for further work is whether this construction is 
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the best choice, or if other schemes (including the symmetric one) are more viable in a numerical 
sense. 

On the other hand, in [52] a complete quantization of the vacuum Bianchi I model was deve- 
loped within the p, scheme. Subsequently, the analysis was also extended to the fi' quantization 
prescription [54]. In this approach, in order to have a concept of evolution, the role of the clock 
must be fulfilled by some geometric degree of freedom. In the context of the quantization given 
in [52], two different variables (a triad coefficients and its conjugate momentum) have been 
tested numerically [53] . This numerical analysis has been restricted to Gaussian states and, due 
to complicate numerical problems, long term evolutions could not be checked. Therefore, further 
numerical analysis are needed to conclude, for instance, that all semi-classical states remain so 
after the bounce. Even though, the bounce picture was proven to be robust under the new 
interpretation of time. 

4.6 Effective models 

The approach given by the effective equations is based on the so-called geometric quantum 
mechanics. In this context, the quantum Hilbert space is considered to be an infinite dimensional 
phase space with the symplectic form given by the imaginary part of the inner product. Taking 
expectation values one can relate both spaces: each operator on the Hilbert space is mapped 
to a function on the phase space. Then, it can be seen that the commutator operation of the 
Hilbert space is directly related to the Poisson brackets of the phase space. Since the equations 
of motion need to be solved only for the classical degrees of freedom (that are interpreted as 
the expectation values of the observables in the quantum theory) with given quantum geometric 
corrections, one deals with differential equations and not with discrete equations. In particular, 
for the case of homogeneous systems, we only need to deal with ordinary differential equations 
since the unique dependence of all objects will be the time variable. This fact makes easier, 
specially from a computational point of view, to analyze the different models effectively rather 
than facing the full quantum evolution equations. 

As far as we know, there is only one analysis [72] in the literature, corresponding to the 
quantization of FRW model, that has systematically derived the effective Hamiltonian and 
equations of motion by computing the expectation value of the quantum Hamiltonian operator 
for Gaussian coherent states and then performing an asymptotic expansion. For the rest of 
the cases, the effective equation program has essentially consisted on replacing the operators 
in the quantum Hamiltonian constraint by their expectation values. This provides a first order 
approach inasmuch as one disregards all state dependent information. But it gives rise to 
a modified classical Hamiltonian that contains corrections by LQG effects, from which the 
effective equations of motion for different variables can be obtained. The effective evolution of 
these variables is expected to mimic the behavior of the quantum system. In fact, in all the cases 
presented in this section, that the quantum discrete equation has been numerically solved, the 
pure quantum results have been compared with the effective ones, and a very precise agreement 
has been obtained. 

Effective descriptions of quantum FRW are used to intuitively understand different aspects of 
the system and generalize the obtained results through a purely quantum analysis. For instance, 
the analytical form of pcrit has been obtained [5] . In addition, the presence of a short phase of 
superinflation immediately after the bounce was first shown for a scalar matter field in [15]. This 
result was afterwards generalized to other matter models [67]. During this superinflationary 
phase both the Hubble parameter and its time derivative are positive; thus the acceleration 
of the universe is faster than in a de-Sitter spacetime. Furthermore, the resolution of strong 
singularities has also been shown throughout this effective analysis [68]. A particularly important 
question is whether LQG may leave any testable footprints in the GMB. This issue has been 
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analyzed also by effective means in several articles (see, for instance, [73, 45, 66, 43, 57, 20]) 
with a wide variety of assumptions, but yet no final consensus has been achieved. 

On the other hand, effective equations have been applied to the analysis of more complicated 
models, for which a complete quantum treatment is still out of reach. Regarding the anisotropic 
models, in [23] it has been shown that the chaotic behavior of the classical Bianchi IX model 
near the singularity is reduced by quantum effects for a /io scheme. Effective equations for this 
model in the p! scheme have already been proposed [75], but a numerical analysis is still missing. 

The effective dynamics of the vacuum Bianchi I model was first analyzed in [41] in the context 
of jiQ scheme. This case is simple enough so that the analytical solution can be achieved, which 
clearly shows the resolution of the classical singularity since the volume is not allowed to vanish. 
For generic matter content, and in the context of the jl scheme, the effective study was carried 
out in [36]. In particular, the vacuum, scalar field and radiation matter contents were analyzed. 
For the particular case of a scalar matter field, in [35] a comparison between the Jl and Ji' schemes 
was presented. This analysis showed that in both cases the classical singularities are resolved 
and replaced by big bounces (one for each spatial direction) but certain key details differ. In 
the scheme these bounces occur when the corresponding directional densities pi := p^/p^, 
which depend on the fiducial cell, reach certain critical value. Whereas, for the fl' case, the 
characterization of the bounce is given by the critical value of the matter density. Thus, in the 
last case, all three bounces take place almost simultaneously. The classical collapsing universe 
that is obtained "on the other side of the bounces" also depends on the quantization scheme. 

There are also several effective analysis of the interior of the Schwarzschild black hole and 
even if, as far as we know, none of them is based on a rigorous quantum theory, appealing 
results have been obtained. For instance, in [13] the singularity resolution was already obtained 
for both Ho and p,' quantization schemes. Constructing on these results, in a manner parallel 
to the above-commented Bianchi I case [35] , Chiou performed a systematic comparison between 
the physical consequences of p, and p' schemes for this model [37]. In both cases the classical 
singularity is resolved but the final picture is quite different. Taking into account an extended 
Schwarzschild spacetime, in the p scheme quantum bounces replace the classical white and black 
hole singularities and work as a bridge between their corresponding regions. On the other hand, 
in the p' case, the classical singularity is resolved, which implies that the even horizon is diffused. 
This produces a baby black hole with a very reduced mass in the consecutive classical cycle. 
A process that continues until the spacetime enters a highly quantum regime. 

Regarding cosmological inhomogeneous models, the vacuum Gowdy spacetime (with a three- 
torus spatial topology) has been numerically analyzed in [29] for the p scheme. This effective 
treatment is based on the hybrid quantization of the mentioned model [50], where different 
degrees of freedom are treated on a distinct footing. The homogeneous degrees of freedom 
are quantized polymerically, whereas a regular Fock quantization is developed for the inhomo- 
geneities. The main result of this analysis is that the bounce picture is robust in the presence 
of inhomogeneities. In addition, a Monte Carlo analysis showed that the amplitude of the inho- 
mogeneities is statistically conserved through the bounce for highly inhomogeneous universes, 
whereas it is amplified for almost homogeneous cases. 

In order to finalize, we would like to briefly mention that a different effective method (the 
so-called truncation method [8]) has been used in quantum cosmology [26, 25]. Following the 
geometric quantization mentioned above, a state can be completely parameterized by the expec- 
tation values of the basic operators, like q and p in quantum mechanics, and an infinite number of 
moments (expectation values of 4""^^ for all a and b) . The equations of motion for all these vari- 
ables, which constitute an infinite system of coupled differential equations, is then completely 
equivalent to the quantum dynamics of the corresponding wave function. Nonetheless, in certain 
(semiclassical) situations, this infinite system can be truncated by keeping terms up to a given 
order (defined as the sum a + b). In this way, it is possible to analyze the back-reaction of the 
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high-order moments on the trajectories of the expectation values. For instance, this approach 
has been used in [24] to obtain and analyze the effective equations (for both WdW and loop 
quantizations) of cosmological homogeneous models with a massive or interacting scalar field 
at second order. In contrast to the free scalar field case, the back reaction is very relevant in 
this model since the classical equations for the expectation values receive corrections due to the 
back-reaction of fluctuations. The LRS Bianchi I model was also analyzed in [27]. In this case, 
the matter content was chosen as an isotropic scalar field with a negative energy density given 
in terms of the anisotropic scale factors. This exotic matter field makes the system treatable in 
terms of physical coherent states. Therefore, a comparison between the effective and the purely 
quantum system was performed and both approaches showed an excellent agreement. Finally, 
in [19], in order to systematize the computation of equations of motion for the moments at very 
high order, powerful algebraic computational tools have been constructed. Furthermore, these 
tools have then been applied to the study the quantum back reaction on a homogeneous universe 
with positive cosmological constant. 

5 Summary 

In this article we reviewed a variety of numerical techniques and their application to several 
relevant LQC models. We commented on the features of the quantum Hamiltonian constraint 
in the context these models that make these techniques readily applicable and extremely useful, 
especially for exploring different quantization schemes and examining the semi-classical limit. 
More specifically, inspired by its role in numerical analysis, we reviewed the concept and appli- 
cation of von Neumann stability to several LQC models. We also reviewed the significant role 
played by lattice refinement in the context of the semi-classical behavior of these models. This 
article also includes a survey of a wide variety of isotropic and anisotropic LQC models and how 
these numerical techniques have contributed towards a deeper understanding of their features, 
and a discussion of open issues that require further work to advance this emerging field. 

A Stability analysis of the closed FRW model 

As an example of stability analysis, we consider here the case of the k = 1 FRW model, with two 
incarnations of the quantization method - the work of Green and Unruh [44] with the old quanti- 
zation scheme, and Ashtekar, Pawlowski, Singh and Vandersloot [7] with the improved quantiza- 
tion method. In the original Green and Unruh paper, the authors found that their Hamiltonian 
constraint generically led to solutions that increased in size without bound for large values of 
the triad parameter p, obviously not in tune with the classical recollapse expected for a closed 
isotropic model. Their discussion of these results notes the absence of an inner product in their 
framework to eliminate such unphysical states, but they posited that the lack of a recollapse indi- 
cated the LQC program has less contact with the full LQG theory than previous results had indi- 
cated. However, the results of Ashtekar et al., using a Hamiltonian constraint with an improved 
quantization scheme with lattice refinement, put to rest those doubts by finding both a bounce 
in the quantum regime near the classical singularity, and a recollapse at the appropriate classi- 
cal scale. In their comparison of the results in the two papers (Section VIIB in [7]), they point 
out two differences between the constraints - the self-adjoint nature of the improved quantum 
constraint, and the use of a scalar field by Ashtekar et al. as a time parameter for the evolution. 

However, we show here that the Green and Unruh constraint suffers from another problem 
beyond this - namely the presence of an unstable "numerical" mode that grows without bound, 
leading to the lack of classical recollapse. We now derive this using von Neumann stability 
analysis, and show this mode does not exist for the improved quantization scheme of Ashtekar 
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et al. For both models, we make the ansatz for the wavefunction 

Here, m is a placeholder for the appropriate variable of the two models (scaled as appropriate), 
while (j) is the value of the scalar field. For the Green-Umuh model, this is /i, which is pro- 
portional to the triad component p; for the improved model, this is the volume v. With both 
equations in this ansatz, the quantum Hamiltonian constraint is of the form 

C+(m)-i/'m+4 + C'°(m)V'm + C"(m)V'm-4 = to'^B{m)ilJm- 

If in addition, we assume V'm = A™/*^ in the asymptotic limit m — )• oo, where A is a constant, 
then this leads to a quadratic equation in A, namely 

C+\^ + {Cl - u;^Bo)X + Cq = 0, 

where Cq , Cq, Bq and Cq are the limits of the appropriate function for m — t- oo. Using the 
standard formula, this gives a solution 



A — i • 

2Co+ 

If |A| > 1, then we have a solution corresponding to this constant which increases in magnitude 
without bound; this can be easily seen, since in our ansatz, 

hm ( ^1 = A. 

Now we look at the particular functions for each model, and find their corresponding values 
of A. 

For the Green-Unruh model, we have that 

C+{m) = [V{m + 5) - V{m + 3)]e-^'^''°, 

C\m) = -[2 - Afil-f\T^ - T)][V{m + 1) - V{m - 1)], 

C-{m) = [V{m - 3) - V{m - 5)]e^'^'"\ 



with the volume function given by 
|m|/Uo7^p 



V{m) 
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the function dj^i{m) coming from the quantum ambiguity in the inverse volume operator [17], 
defined as 
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j-(j + l)(2i + 1)7^2/ 




(m//o + 2r) 



3/(2-20 



and r = 1/2 indicates the closed model. For large values of m, then the inverse volume op- 
erator matches the classical expression closely, so that dj^i{m) ~ (m/Uo7^p/6)~^/2, Thus, the 
asymptotic expressions for the coefficient functions are given by 
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Note that the scalar field eigenvalue ui drops out, because the inverse volume dependence B{m) ~ 
jg ^Q-^gj- order than that of C^(m) ~ m^/^. This gives values to leading order of 

Thus, one eigenvalue is larger than one in magnitude; using the constant value //q = \/3/4 and 
7 = 0.2735, then |A+| ~ 1.126 while |A_| ~ 0.8884. 

Turning now to the improved constraint of Ashtekar et al., we have 

C+{v) = ^|t; + 2|||?; + l| -\v + 3\\, C- {v) = C+iv-i), 
8 

C^{v) = -C+{v) -C {v) + — l3K SIV? ( - ^ - i-O - j 

X ||w + 1| — |u — 1||, 
B{v)=(^^\\v\\\v + l\^/^ -\v-l\^/^\' 

and fi{v) = 2,^/2/2{K/vf/^. Using the same methodology as before, we have that for large 
magnitude of f , 

= C„- = 5^ + 0(1), CS - .-B, = + 0(m'/3). 

Plugging these into the equation for X± gives A± = 1, and the equation is stable. 
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